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EXECUTIVE  SUMMARY 


OBJECTIVE 

Given  any  linear  system  of  equations  Mx  =  b  in  which  M  is  a  large  sparse 
symmetric  matrix,  provide  a  fully  automated  computer  program  for  generating 
computational  tasks  that  can  be  processed  independently  of  each  other  by  the 
different  processors  of  a  parallel  architecture  computer.  Such  an  automated 
parallelization  tool  is  essential  for  the  effective  applications  of  parallel  architecture 
computers. 

RESULTS 

VVs  have  presented  a  detailed  computer  implementation  of  a  combinatorial 
algorithm  developed  in  part  1  (Kevorkian,  1993)  for  decomposing  a  large  sparse 
symmetric  system  of  equations  Mx  =  b  into  independently  solvable  smaller  tasks 
that  can  be  executed  in  parallel  on  different  processors  of  a  parallel  architecture 
computer.  Also,  we  have  presented  a  procedure  that  uses  the  output  of  the 
computer  program  to  generate  a  block  bordered  diagonal  form  of  matrix  M  that  is 
well-suited  for  sparse  block  factorizations. 
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1.  INTRODUCTION 


The  solution  of  large  sparse  linear  symmetric  systems  of  equations  of  the  type 


Mx  =  b 

forms  the  most  compute-intensive  part  in  several  of  the  Navy's  fundamental  grand 
challenge  problems.  These  problems  include  ocean  basin  scale  modeling,  three- 
dimensional  modeling  of  ocean  acoustic  propagation,  fluid  flow  simulations,  tur¬ 
bulent  combustion,  and  structural  design  of  navy  vehicles,  as  well  as  linear  and 
nonlinear  constrained  optimization  problems  known  as  linear  and  nonlinear  pro¬ 
gramming.  In  most  of  these  applications,  there  is  significant  parallelism  hidden  in 
the  structure  of  the  original  problem.  Therefore,  any  progress  toward  the  efficient 
solution  of  a  general  grand  challenge  problem  on  a  parallel  architecture  computer 
will  require  advanced  computational  tools  that  can  exploit  all  hidden  parallelism. 

In  this  work  we  give  a  complete  computer  implementation  of  an  algorithmic  tooi 
developed  recently  by  Kevorkian  (1993)  for  exploiting  parallelism  hidden  in  the 
sparsity  structure  of  large  sparse  symmetric  matrices  with  regular  and  irregular 
structures.  The  application  of  this  automated  parallelization  tool  to  a  general  sparse 
symmetric  system  of  equations  decomposes  the  original  problem  into  inde¬ 
pendently  solvable  smaller  tasks  for  execution  on  different  processors  of  a  parallel 
architecture  computer. 

This  report  is  organized  as  follows.  In  section  2,  we  present  a  computer  imple¬ 
mentation  of  the  parallelization  tool,  called  "roadmap"  (Kevorkian,  1993),  using  the 
linear  algebra  package  Matlab  (MathWorks,  1990).  In  section  3,  we  present  a 
procedure  that  uses  the  program  roadmap  to  generate  a  permutation  matrix  P  such 
that  PMpT  is  a  block  bordered  diagonal  matrix  satisfying  all  three  properties  of  the 
vertex  partition  computed  in  roadmap,  in  section  4,  we  discuss  future  works 
pertaining  to  experimental  results  obtained  from  the  application  of  roadmap  to  a 
collection  of  test  problems.  Also,  we  discuss  the  development  of  a  recursive  version 
of  roadmap  that  can  exploit  parallelism  not  only  in  the  original  problem,  but  also  in 
all  subsequent  parts  (Schur  complements)  that  usually  result  from  the  solution  of  a 
sparse  system  of  equations  by  Gaussian  elimination. 

2.  IMPLEMENTATION  OF  PARALLELIZATION  TOOL  ROADMAP 

In  this  section  we  present  a  complete  computer  implementation  of  roadmap  by 
using  the  widely  available  linear  algebra  package  Matlab  (MathWorks,  1990). 

The  input  to  roadmap  is  an  n-by-n  structurally  symmetric  sparse  matrix  M  =  [mjj] 
with  the  data  structure  described  as  follows. 

Suppose  G  =  (V,  E)  denotes  the  undirected  graph  of  the  n-by-n  matrix  M.  Then 
the  set  V  consists  of  n  vertices  vi  through  Vn,  with  Vj  representing  row  i  of  matrix  M. 
For  convenience,  let  Vj  =  i,  for  i  =  1  ,...,n.  Then  we  get  adjQVi  =  { j  |  mjj  0  for  all  j  i }, 
and  so  the  array  ADJ(i)  consists  of  the  set  of  column  indices  of  the  nonzero  off- 
diagonal  entries  in  row  i  of  matrix  M.  Thus,  the  problem  of  representing  the  matrix  M 
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by  a  sparse  data  structure  is  one  of  finding  a  compact  way  of  storing  and 
referencing  the  n  arrays  ADJ(1)  through  ADJ(n).  The  sparse  data  structure  we  use 
in  roadmap  was  first  introduced,  implemented,  and  applied  by  Gustavson  (1973). 

Suppose  we  combine  the  n  arrays  ADJ(1)  through  ADJ(n)  to  form  a  single  array 
LADJ  defined  by 


LADJ  =  [ADJ(1),  ADJ(2),  ...  ,  ADJ(n)  ]. 

Also,  let  lADJ  be  an  n+1  single  array 

lADJ  =  [IADJ(1),  1ADJ(2), ... ,  lADJ(n+1)] 
with  elements  defined  by 

i-1 

IADJ{i)=1+ DEGG) ,  i  =  1 . n+1. 

j=i 

Since  DEG(i)  =  lADJ(i)|,  the  integers  lADJ(i)  and  IADJ(i+1)  -1  are  pointers  to  the  first 
and  last  vertices  of  ADJ(i)  on  array  LADJ,  respectively,  and  so  using  standard 
Matlab  notation  (MathWorks,  1990)  we  obtain  the  following  two  equalities; 

ADJ(i)  =  LADJ(IADJ(i):IADJ(i+1)  - 1) , 
and 

DEG(i)  =  IADJ(i+1)  -  lADJ(i),  i  =  1 , ... ,  n. 

Thus  the  data  structure  for  the  zero-nonzero  structure  of  matrix  M  is  a  pair  of  single 
arrays  (lADJ,  LADJ),  where  lADJ  consists  of  n+1  pointers  and  LADJ  consists  of  2|E| 
column  indices.  George  and  Liu  (1981)  refer  to  the  pair  of  arrays  (lADJ,  LADJ)  as 
the  adjacency  structure  pair  of  matrix  M. 

The  entire  roadmap  program  is  given  in  figures  1  through  8.  In  writing  this  pro¬ 
gram,  our  main  objective  was  to  establish  a  concise  one-to-one  correspondence 
between  the  implementation  and  the  high-level  language  used  in  Kevorkian  (1993) 
so  that  all  theoretical  parts  of  the  work  can  be  conveniently  traced  and  studied.  All 
parameters,  variables  and  arrays  used  in  the  roadmap  are  described  in  detail  in 
Kevorkian  (1993). 

The  procedure  "sprsdata"  called  in  roadmap  queries  the  user  for  the  dimension 
n  of  matrix  M  and  requires  from  the  user  the  adjacency  structure  pair  (lADJ,  LADJ) 
of  M.  Procedure  sprsdata  also  computes  the  single  array  DEG  and  monitors  the 
integrity  of  the  input  data  by  checking  the  correctness  of  three  relations.  These  are 
as  follows: 


llADJ|  =  n+1, 
|LADJ|  =  IADJ(|IADJ|)-1, 
LADJ(i)  <n. 


foralli<|LADJ|. 
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If  any  of  these  three  relations  is  violated,  the  program  halts  and  prints  the  message 
"check  data."  The  program  also  halts  if  the  following  equality  holds 

|LADJ|  =  nx(n-1) 

since  G  is  a  clique  in  this  a  case;  and  so  no  sparsity  will  exist  in  G. 

Procedure  "search"  called  in  roadmap  computes  the  set  of  vertices  S.  If  the  set 
S  is  nonempty,  then  roadmap  calls  procedure  "dfs"  to  compute  the  connected  com¬ 
ponents  of  the  induced  subgraph  G(V-S).  Othenwise,  roadmap  calls  procedure 
"cliques"  to  compute  independent  cliques  in  the  original  graph  G. 

The  procedure  "classify"  called  in  roadmap  categorizes  the  connected  compo¬ 
nents  of  G(V-S)  into  cliques  and  noncliques  using  parameters  computed  in  dfs.  If  a 
connected  component  G(U)  is  a  clique,  then  procedure  classify  uses  results  estab¬ 
lished  in  Kevorkian  (1993)  to  classify  G{U)  into  one  of  four  distinct  types  of  cliques. 

If  G(U)  is  not  a  clique,  then  procedure  classify  calls  procedure  cliques  to  compute 
independent  cliques  in  the  nonclique  connected  component  G(U). 

The  clique  connected  components  of  induced  subgraph  G(\/-S)  combined  with 
the  independent  cliques  computed  in  each  of  the  nonclique  connected  compo¬ 
nents  of  G(V-S)  form  the  independently  solvable  smaller  tasks  that  can  be  exe¬ 
cuted  in  parallel  on  different  processors  of  a  parallel  architecture  computer. 

3.  BLOCK  BORDERED  DIAGONAL  FORMS  USING  ROADMAP 

Procedure  bbdf  given  in  figure  9  uses  the  program  roadmap  to  generate  a 
permutation  matrix  P  such  that  PMP^  is  a  block  bordered  diagonal  matrix  satisfying 

all  three  properties  of  the  vertex  partition  n*  =  {Vi,  V2 . Vr,  S*).  These  properties 

are  briefly  stated  in  procedure  roadmap  shown  in  figure  1 ,  and  covered  in  more 
depth  in  Kevorkian  (1993).  By  the  first  two  properties  of  the  vertex  partition,  PMPl’  is 
an  (r+1)-by-(r+1)  block  bordered  diagonal  matrix  such  that  each  of  the  leading  r 
diagonal  blocks  is  a  full  matrix.  By  the  third  property,  every  principal  submatrix  in  M 
that  corresponds  to  an  interior  clique  in  the  graph  of  M  is  a  diagonal  block  in  PMPT 
In  part  1  of  this  work  (Kevorkian,  1993),  we  have  shown  that  the  symbolic 
factorization  of  a  matrix  corresponding  to  an  interior  clique  does  not  produce  any 
fill-in.  Block  bordered  diagonal  forms  computed  by  the  program  roadmap  are  thus 
well-suited  for  sparse  block  factorizations. 

The  procedure  bbdf  consists  of  three  distinct  parts.  The  first  part  uses  the 
adjacency  structure  pair  (lADJ,  LADJ)  to  generate  a  Boolean  form  of  matrix  M.  The 
second  part  uses  the  array  QUEUE  and  properties  of  arrays  SN  and  VC  to  compute 
the  vertex  ordering  a.  Permuting  the  rows  of  M  in  the  sequence  given  by  a  produces 
the  block  bordered  diagonal  form  of  matrix  M.  The  third  and  last  part  of  procedure 
bbdf  uses  the  signs  of  the  elements  placed  on  IQUEUE  to  generate  an  r+1  array 
such  that  the  ith  element  is  a  pointer  to  the  first  row  of  the  ith  diagonal  block  in 
PMP''’.  While  modifying  the  array  IQUEUE  in  procedure  bbdf,  we  make  sure  that  the 
pointers  to  the  starting  vertices  of  independent  cliques  computed  in  cliques  retain 
their  original  negative  signs.  This  way  we  are  able  to  use  the  array  TYPE  computed 
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in  roadmap  to  relate  each  diagonal  block  in  PMP''"  to  the  type  of  clique  it  is 
associated  with  in  G. 

The  first  two  parts  of  procedure  bbdf  are  straightfonward  and  easy  to  follow.  The 
third  part  is  more  complicated,  requiring  the  following  property  of  array  IQUEUE  for 
correctness. 

Lemma  1.  Let  S  be  the  set  of  vertices  computed  in  procedure  search.  Let  i  and  j 
be  any  two  consecutive  elements  on  the  array  IQUEUE  at  the  completion  of 
roadmap.  Then  the  set  of  vertices  on  QUEUE(li|:|jl-1)  is  a  subset  of  the  vertex  set 
S*-  S  if  and  only  if  i  <  0  and  j  >  0. 

Proof.  Let  i  and  j  be  any  two  consecutive  elements  on  array  IQUEUE.  Then  one 
of  the  following  two  cases  must  hold. 

Case  f.  i  >  0.  Then  by  the  construction  of  procedure  dfs,  there  is  a  connected 
component  G(U)  of  G(V-S)  such  that  the  vertex  v  =  QUEUE(i)  is  the  starting  vertex 
of  the  connected  component  G(U)  computed  in  dfs.  Suppose  j  >  0.  Then  by  the  last 
statement  in  procedure  maxclq,  it  follows  that  no  call  was  made  to  procedure 
cliques  in  classify  at  the  completion  of  the  connected  component  G(U).  This  means 
that  the  connected  component  G(U)  is  a  clique,  and  so  no  part  of  the  set  of  vertices 
on  QUEUE(i;j-1)  is  in  the  set  S**S.  Now  suppose  j  <  0.  Since  i  >  0,  the  integer  j 
must  be  the  first  element  added  to  array  IQUEUE  in  procedure  maxclq  at  the 
completion  of  connected  component  G(U).  Thus  G(U)  is  not  a  clique  and, 
furthermore,  the  vertices  on  QUEUE(l:|j|-1)  comprise  the  vertex  set  in  the  first 
independent  clique  computed  by  procedure  cliques  in  G(U).  As  a  result,  no  part  of 
the  set  of  vertices  on  QUEUE(i:lj|-1)  can  be  in  the  set  S*-S.  Hence  for  any  i  >  0  and 
j  >  0  or  any  i  >  0  and  j  <  0,  no  part  of  the  set  of  vertices  on  QUEUE(i:lj|-1 )  is  in  the 
set  S**S. 

Case  2.  i  <  0.  Assume  j  <  0.  Then  by  the  construction  of  procedure  cliques,  the 
vertex  v  =  QUEUE{|i|)  is  the  starting  vertex  of  an  independent  clique  computed  by 
cliques  in  some  nonclique  connected  component  G(U)  of  G(V-S),  whereas  the 
vertex  w  =  QUEUE(lj|)  is  either  the  staring  vertex  of  another  independent  clique  in 
G(U)  or  is  the  starting  vertex  of  a  separator  computed  by  cliques  in  G(U).  Therefore 
for  the  case  where  j  <  0,  the  vertices  on  QUEUE(|i|:Ij|-1 )  comprise  the  vertex  set  in 
an  independent  clique  computed  by  procedure  cliques  in  G(U)  and  so  no  part  of 
the  set  of  vertices  on  QUEUE(|i|:|j|-1)  can  be  in  the  set  S*-S.  Finally,  suppose  j  >  0. 
Then  v  is  the  starting  vertex  of  the  separator  computed  by  cliques  in  G(U),  which 
means  that  the  set  of  vertices  on  QUEUE{|il;j-1)  is  in  the  set  S*-S. 

This  completes  the  proof. 
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%  procedure  roadmap 
% 

****************************************************************************************** 

%  Given  a  sparse  symmetric  matrix  M,  roadmap  computes  in  the  undirected  graph  * 
%  G  =  (V,  E)  of  M  a  vertex  partition  n*  =  (Vi,  V2, ...  ,  Vr,  S*)  satisfying  the  following  * 
%  three  properties: 

%  (a)  for  any  two  distinct  elements  V;  and  Vj  of  the  partition,  no  vertex  in  Vj  is 

%  adjacent  to  a  vertex  in  Vj, 

%  (b)  every  element  Vj  of  the  partition  induces  a  clique, 

%  (c)  interior  of  every  clique  in  G  is  an  element  of  the  partition. 

%  Program  roadmap  computes  the  vertex  partition  n*  in  linear  time. 

A**********-****^************************************************#***********#****#******** 


% 

global  lADJ  LADJ  LEAF  NGU  RANKE  SN  TEST  U  VC 


sprsdata 
QUEUE  =  []; 
IQUEUE  =  D: 

TYPE  =  D; 

LEAF  =  0; 

VC  =  zeros(1,n); 
SN  =  zeros(1,n); 
TEST  =  zeros(1,n); 
search 

if  any(SN  ==  1) 
dfs 


%  user  provided  input  (n,  lADJ  and  LADJ) 
%  store  connected  components  of  G(V-S) 
%  pointers  to  "roots"  and  "end"  of  QUEUE 
%  classify  "type"  of  connected  component 
%  pointer  to  last  vertex  placed  on  QUEUE 
%  mark  all  vertices  "new" 

%  initialize  separator  numbers  to  zero 
%  initialize  Boolean  array  TEST  to  zero 
%  compute  the  set  S  =  { v  |  SN{v)  =  1 } 

%  if  S  is  nonempty  then 
%  compute  connected  components 


else 

%  else 

R00T  =  1: 

% 

pointer  to  root  vertex  of  V 

LEAF  =  n; 

% 

pointer  to  end  vertex  of  V 

QUEUE  =  [ROOTiLEAF]: 

% 

place  entire  vertex  set  V  on  QUEUE 

IQUEUE  =  [ROOT]: 

% 

pointer  to  root  vertex  on  QUEUE 

TYPE  =  [Oj; 

% 

G  is  a  regular  graph  (&  not  a  clique) 

cliques 

% 

compute  independent  cliques  of  G 

end 

% 

end 

IQUEUE  =  [IQUEUE,  LEAF+tj; 

%  pointer  to  end  of  QUEUE 

Figure  1.  Procedure  roadmap. 
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%  procedure  sprsdata 
% 

^******«  «***«  «  ************************************************************************* 
%  This  p.  .^cedure  performs  the  following  four  tasks; 

%  (a)  queries  the  user  for  size  of  matrix  fn); 

%  (b)  requires  from  user  the  adjacency  structure  pair  (lADJ,  LADJ)  of  matrix;  * 

%  (c)  tests  correctness  of  input  data; 

%  (d)  computes  the  single  array  DEG  (degrees  of  vertices). 

%************* . * . * . ***** . ****** . ***** . * . 

% 

%  query  user  for  size  of  matrix  (n) 
n  =  input('Enter  n  (if  n  <=  1 ,  program  quits): '); 
if  n  <=  1 ,  break,  end 
end 

%  provide  adjacency  structure  pair  (lADJ,  LADJ) 
if  n  ==  21 

%  Illustrative  example  used  in  Kevorkian  (1993) 

lADJ  =  [1  6  8  1 1  14  23  27  32  41  44  47  5l  55  61  67  70  75  82  86  90  93  97); 
LADJ  =  [578  16  17  13  17  4  1420  3  1320  1  678  14  16  17  18  19]; 

LADJ  =  [LADJ,  5  14  18  19  1  5  8  16  17  1  5  7  11  12  13  16  17  21  10  15  20]; 
LADJ  =  [LADJ,  9  14  15  8  12  13  21  8  11  17  21  2  48  11  15  21]; 

LADJ  =  [LADJ,  3  56  10  18  19  9  1013  1  5  78  17]; 

LADJ  =  [LADJ,  1  2  5  78  12  16  5  6  14  19  5  6  14  18  349  8  11  12  13]; 
end 

%  test  correctness  of  input  data 
RANKI  =  length(IADJ); 

RANKL  =  length(LADJ); 
if  RANKI  ~=  n+1 

dispC  lADJ  does  not  have  correct  length.  Check  data.'); 
break,  end 
end 

if  RANKL  ~=  IADJ(RANKI)-1 

dispC  LADJ  does  not  have  correct  length.  Check  data.’); 
break,  end 
end 

if  any(LADJ  >  n) 

dispC  LADJ  contains  a  column  index  >  n.  Check  data.'); 
break,  end 
end 

if  RANKL  ==  n*(n-1) 

dispC  Off-diagonal  part  of  M  is  full  and  so  G  is  a  clique.'); 
break,  end 
end 

%  compute  degrees  of  vertices 
DEG=  IADJ(2:n+1)-IADJ(1  :n); 


Figure  2.  Procedure  sprsdata. 
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%  procedure  search 
% 

0/  ************************************************************************************* 

/o 

%  This  procedure  computes  in  G  =  (V,  E)  a  set  of  vertices  S  defined  by 
% 

%  S  =  { V  I  there  exists  in  E  an  edge  (v,  w)  with  DEG(v)  >  DEG(w)} . 

% 

%  The  set  S  has  the  property  that  every  interior  clique  in  G  is  a  connected 
%  component  of  induced  subgraph  G(V-S):  Corollary  3.1  in  Kevorkian  (1993).  * 
(^********«*********************** ******************************************************** 
% 

for  V  =  1  :n-1  %  for  each  vertex  v  in  V  -  (n)  do 


vc(v)  =  i: 

% 

mark  vertex  v  "old" 

forw  =  LADJ{IADJ{v):IADJ(v+1)-1) 

% 

for  ail  w  adjacent  to  v  do 

if  VC(w)  ==  0 

% 

if  w  is  marked  "new"  then 

if  DEG(v)  ~=  DEG(w) 

% 

if  DEG(v)  ^  DEG(w)  then 

if  DEG(v)  <  DEG(w) 

% 

if  DEG(v)  <  DEG(w)  then 

SN(w)  =  1 : 

% 

w  is  in  S 

else 

% 

else 

SN(v)  =  1; 

% 

V  is  in  S 

end 

% 

end 

end 

% 

end 

end 

% 

end 

end 

% 

end 

end  %  end 

VC  =  zeros(1  ,n);  %  mark  all  vertices  "new" 


Figure  3.  Procedure  search. 
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%  procedure  dfs 
% 

(^«**4************************************************************************************** 

%  This  procedure  computes  all  connected  components  of  induced  subgraph 
%  G(V-S).  We  use  G(U)  to  denote  a  connected  component  computed  in  dfs. 

{^********** ************************************************************************ ******** 
% 


for  V  =  1  :n 
if  SN(v)  ==  0 
if  VC(v)  ==  0 
VC(v)  =  1 

QUEUE  =  [QUEUE, V]: 

LEAF  =  LEAF+1: 

ROOT  =  LEAF; 

IQUEUE  =  [IQUEUE,  ROOT]; 
RANKE  =  0; 

NGU  =  []; 
cmponent(v) 

RANKU  =  LEAF-ROOT+1; 
RANKN  « length(NGU); 
TEST(NGU)  =  zeros(1 , RANKN); 
classify 
end 
end 
end 


%  for  all  V  in  V  do 
%  if  V  is  in  V/-S  then 
%  if  V  is  marked  "new"  then 

%  mark  v  "old" 

%  V  is  the  root  vertex  of  G(U) 

%  pointer  to  end  vertex 

%  pointer  to  root  vertex 

%  add  root  pointer  to  IQUEUE 

%  count  edges  visited  in  dfs 

%  neighborhood  of  U  in  G 

%  compute  G(U) 

%  compute  size  of  U 

%  compute  size  of  NGU 

%  update  Boolean  array  TEST 

%  identify  "type"  of  component 

%  end 

%  end 

%  end 


Figure  4.  Procedure  dfs. 
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%  function  cmponent(v) 


(^** *************************************************************************** *********** 


%  This  function  computes  neighborhood  NGU  of  U  while  computing  G(U) 

O/  **************************************************************************************** 

/o 


function  cmponent(v) 
for  w  =  LADJ(IADJ(v);IADJ(v+1  )-1 ) 
if  SN(w)  ==  0 

RANKE  =  RANKE+1 ; 
if  VC(w)  ==  0 
VC(w)  =  1: 

QUEUE  =  [QUEUE,  w]  ; 
LEAF  =  LEAF+1 ; 
cmponent(w) 
end 

^'^?f  TEST(w)  ==  0 

NGU  =  [NGU,  w]: 
TEST(w)  =  1  ; 
end 
end 
end 


%  declare  cmponent(v)  as  function  file 
%  for  all  w  adjacent  to  v  do 
%  if  w  is  in  V-S  then 

%  account  for  visited  edge  (v,w) 

%  if  w  is  marked  "new"  then 

%  mark  w  "old" 

%  add  w  to  QUEUE 

%  update  pointer  to  end  vertex 

%  call  cmponent(w) 

%  end 

%  else 

%  if  w  is  not  on  NGU  then 

%  add  w  to  NGU 

%  mark  w  as  vertex  on  NGU 

%  end 

%  end 

%  end 


Figure  5.  Function  cmponent. 
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%  procedure  classify 
% 

0/  ****************************************************************************************** 

/o 

%  This  procedure  categorizes  the  connected  components  of  G(V-S)  into  cliques 
%  and  noncliques  using  the  algebraic  relation  RANKE  =  RANKU*(RANKU'1).  * 
%  Subsequently,  all  clique  connected  components  in  G(V-S)  are  classified  into 
%  the  four  types  of  cliques  Ci  through  C4  using  Corollaries  4.2  and  4.3  in 
%  Kevorkian  (1993). 

****************************************************************************************** 


/rANKE  ==  RANKU*(RANKU-1 )  %  if  G(U)  is  a  clique  then 


if  RANKN  ==  1 

% 

if  Cor.  4.2  holds  then 

TYPE  « [TYPE,  1]; 

% 

G(U)  is  an  interior  clique 

else 

% 

else 

R  =  RANKN+RANKU-1; 

% 

compute  integer  R 

If  any(DEG(QUEUE(ROOT:LEAF))  R) 

% 

If  Cor.  4.3  does  not  hold 

TYPE  =  [TYPE,  3]; 

% 

G(U)  is  not  an  si  clique 

else 

% 

else 

TYPE  =  [TYPE,  -2]: 

% 

G(U)  is  an  si  clique 

end 

% 

end 

end 

% 

end 

else 

%  else 

TYPE  =  [TYPE,  0]: 

% 

G(U)  is  not  a  clique 

VC(QUEUE(ROOT:LEAF))=zeros(1,RANKU): 

% 

mark  vertices  in  G(U)  "new" 

cliques 

% 

compute  ind.  cliques  in  G(U) 

end  %  end 


Figure  6.  Procedure  classify. 


%  procedure  cliques 
% 

<^** ***************************************************************************** ********* 

%  This  procedure  computes  in  G(U)  independent  cliques  G(Ui),  G(U2) . and  * 

%  their  neighborhoods  N(Ui),  N(U2) .  in  G(U)  such  that 

%  G(Ui)  is  maximal  in  G(U), 

%  G(U2)  is  maximal  in  G(U  -  Ui  -  N(Ui)), 

%  G(U3)  is  maximal  in  G(U  -  Ui  -  U2  -  N(Ui)  -  N(U2)), 

%  and  so  forth.  The  sets  Ui ,  U2, ... ,  are  placed  on  array  CLQS  in  that  order,  * 
%  while  all  neighborhoods  are  placed  on  array  NBRS.  At  the  completion  of 
%  cliques  the  vertex  set  U  on  QUEUE  is  replaced  by  the  array  [CLQS, NBRS]. 

<^***ii-**************************************************  *******************  *************** 
% 

CLQS  =  []:  %  stores  independent  cliques 

NBRS  =  0:  %  stores  neighborhoods  of  ind.  cliques 

CLQROOT  =  ROOT ;  %  pointers  to  ind.  cliques  on  CLQS 

TAIL  =  0;  %  pointer  to  last  vertex  on  CLQS 

for  V  =  QUEUE(ROOT :LEAF)  %  for  each  vertex  v  in  G(U)  do 


if  VC(v)  ==  0 

% 

if  V  is  marked  "new"  then 

vc(v)  =  i: 

% 

mark  v  "old" 

ADJCNT  =  n: 

% 

adj(v)  in  G(U-CLQS-NBRS) 

maxclq 

% 

compute  a  maximal  clique 

end 

% 

end 

end 

%  end 

QUEUE(ROOT:LEAF)  =  [CLQS,  NBRS);  %  replace  U  on  QUEUE  by  [CLQS, NBRS] 

Figure  7.  Procedure  cliques. 
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%  procedure  maxclq 
% 

<^*********************** ********************************************************* *********** 

%  This  procedure  computes  in  induced  subgraph  G(U-CLQS-NBRS)  a  maximal 

%  clique  G(C)  with  starting  vertex  v  (occasionally  called  the  root  vertex). 
******************************************************************************************* 


% 


for  w  =  LADJ(IADJ(v):IADJ(v+1  )-1 ) 

%  for  all  w  adjacent  to  v  do 

if  SN(w)  ==  0 

% 

if  w  is  on  U  but  not  on  NBRS  then 

if  VC(w)  ==  0 

% 

if  w  is  marked  "new"  then 

ADJCNT  =  [ADJCNT,  w]; 

% 

add  w  to  ADJCNT 

else 

% 

else 

NBRS  =  [NBRS,v]: 

% 

reject  v  as  starting  vertex 

SN(v)  =  1; 

% 

avoid  duplicates  of  v 

return 

% 

return  to  cliques 

end 

% 

end 

end 

% 

end 

end 

%  end 

CLQS  =  [CLQS,  V]: 

%C 

=  [v]:  (v  is  starting  vertex  of  G(C)) 

TEST(v)  =  1; 

%  set  TEST(v)  =  1 

RANKC  =  1  ; 

%  RANKC  =  |C| 

for  u  =  ADJCNT 

%  for  each  vertex  u  on  ADJCNT  do 

VC(u)  =  1: 

% 

mark  u  "old" 

w  =  LADJ(IADJ(u):IADJ(u+1)-1); 

% 

w  =  ADJ(u) 

COUNT  =  length(find(TEST(w)  ==  1)); 

% 

COUNT  =  |ADJ(u)  n  C| 

if  COUNT  ==  RANKC 

% 

if  u  is  adjacent  to  all  w  in  C  then 

CLQS  =  [CLQS,  u]; 

% 

C  =  [C,  u] 

TEST(u)=  1; 

% 

set  TEST(u)  =  1 

RANKC  =  RANKC+1 ; 

% 

update  size  of  C 

else 

% 

else 

NBRS  =  [NBRS,  u]; 

% 

u  is  in  neighborhood  of  C 

SN(u)  =  1 ; 

% 

avoid  duplicates  of  u 

end 

% 

end 

end  %  end 

HEAD  =  TAIL+1 ;  %  pointer  to  C(1 )  on  CLQS 

TAIL  =  TAIL+RANKC;  %  pointer  to  C(lCl)  on  CLQS 

TEST(CLQS(HEAD:TAIL))=zeros(1  ,RANKC):  %  reset  TEST  to  zero 
CLQROOT  =  CLQROOT+RANKC;  %  pointer  to  next  starting  vertex 

IQUEUE  =  [IQUEUE,  -CLQROOT];  %  add  pointer  (negated)  to  IQUEUE 


Figure  8.  Procedure  maxclq. 
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%  procedure  bbdf 
% 

<^************ ******************************************************* **«*««***************i^* 

%  This  procedure  uses  the  output  of  program  roadmap  to  generate  a  permutation  * 
%  matrix  P  such  that  RMP"*"  is  a  block  bordered  diagonal  matrix  satisfying  all  three  * 
%  properties  of  vertex  partition  n*  =  (Vi,  V2, ... ,  Vr,  S*). 
<;^************************«*************«******************************«******************** 
% 


M  =  zeros(n,n): 
for  V  =  1  :n 
M(v,v)  =  1 ; 

w  =  LADJ(IADJ(v):IADJ(v+1)-1); 

M(v,w)  =  1 ; 
end 
% 

S  =  QUEUE(find(SN(QUEUE)  ==  1); 

QUEUE  =  QUEUE(find(SN(QUEUE)  ==0)); 
QUEUE  =  (QUEUE, S.  find(SN==1  &  VC==0)] 
M  =  M{QUEUE,  QUEUE); 

% 

if  length(S)  >  0 
SUM  =  0; 
k  =  2' 

i  =  IQUEUE(2): 

for  j  =  IQUEUE{3:length{IQUEUE)) 
if  i  <  0  &  j  >  0 

SUM  =  SUM+i+j: 
else 

IQUEUE(k)  =  i-sign(i)*SUM; 
k  =  k+1 : 
end 
i=j: 
end 
end 


%  use  (lADJ,  LADJ)  to  compute  M 
%  create  an  n  by  n  zero  matrix  M 
%  for  all  V  in  V  do 
%  set  M(v,v)  =  1 
%  w  =  vertices  adjacent  to  v 
%  set  M(v,w)  =  1 
%  end 

%  use  array  QUEUE  to  compute  P 
%  compute  bordered  part  on  QUEUE 
%  compute  block  diagonal  part 
%  vertex  ordering  a  (placed  on  QUEUE) 
%  block  bordered  diagonal  form  of  M 
%  modify  array  IQUEUE 
%  if  |S|  >  0  then 

%  sum  sizes  of  bordered  parts 
%  pointer  for  elements  on  IQUEUE 

%  second  element  on  IQUEUE 

%  for  j  =  IQUEUE(3:|IQUEUE|)  do 
%  if  i  <  0  and  j  >  0  then 
%  update  SUM 

%  ^'^update  IQUEUE 

%  increment  k 

%  end 

%  i  =  jth  element  on  IQUEUE 

%  end 
%  end 


Figure  9.  Procedure  bbdf. 
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4.  FUTURE  WORKS 


Experimental  results  obtained  from  the  application  of  roadmap  to  a  standard  set 
of  test  problems  including  the  Harwell-Boeing  collection  and  a  set  of  matrices 
arising  from  linear  and  nonlinear  programming  optimization  problems  will  be 
reported  in  Kevorkian  (in  preparation-b). 

Also,  we  are  currently  working  on  an  extension  of  roadmap  (Kevorkian,  in 
preparation-a)  that  exploits  parallelism  in  the  original  problem  as  well  as  sub¬ 
sequent  Schur  complements  until  no  further  parallelism  remains  to  exploit.  Such  a 
recursive  version  of  roadmap  will  be  ideally  suited  for  problems  in  which  the 
original  matrix  and  the  Schur  complement  of  the  pivot  block  selected  by  roadmap 
are  very  large  sparse  matrices.  Large  sparse  problems  with  large  Schur  com¬ 
plements  are  frequently  encountered  in  linear  and  nonlinear  programming  opti¬ 
mization  problems  (Kevorkian,  1993). 

5.  CONCLUSION 

We  have  presented  a  detailed  computer  implementation  of  a  linear-time  paral¬ 
lelization  tool  that  automatically  decomposes  a  large  arbitrary  sparse  symmetric 
system  of  equations  into  independently  solvable  smaller  tasks  for  execution  on 
different  processors  of  a  parallel  architecture  computer. 
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